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Abstract.  Three-dimensional  simulation  of  detonation  shock  motion  in  a  condensed 
explosive,  modeled  with  a  reactive  flow  in  a  moderately  complex  geometry  can  require 
enormous  amounts  of  computer  time,  since  the  reaction  zone  behind  the  leading  shock  is 
extremely  thin  compared  to  the  overall  dimensions  of  the  computational  domain.  Therefore 
algorithms  such  as  program  burn,  pre-compute  the  detonation  shock  arrival  time,  and  then 
essentially  use  a  delta  function  model  to  release  the  heat  of  detonation  into  few 
computations  cells  near  the  shock  wave.  Previous  validation  efforts  that  use  a  program  bum 
algorithm  based  on  detonation  shock  dynamics  (DSD)  highlighted  the  prescription  of  the 
initial  detonation  shock  for  the  simulation,  in  a  manner  that  is  self-consistent  with  the 
theory  and  physical  experiment.  In  this  paper  we  use  a  combination  of  theory,  direct  multi¬ 
material  simulation  and  experiment  to  determine  a  critical  radius  for  the  starting  detonation 
shock  shape. 


INTRODUCTION 

In  a  recent  paper  [1],  a  “passover”  experiment 
was  designed  with  a  complex  detonation 
transient.  A  charge  with  an  embedded  lead  disk  in 
a  right  circular  cylindrical  charge  of  HMX-based, 
PBX-9501  was  initiated  from  the  bottom.  A 
detailed  comparison  of  the  motion  of  the 
detonation  shock  was  made  between  experiment, 
prediction  of  the  asymptotic  theory  of  detonation 
shock  dynamic  (DSD-theory)  and  a  multi¬ 
material  reactive  flow  simulation  as  described  in 
[2].  A  new  wide-ranging  equation  of  state  (EOS) 
and  rate  law  [3]  was  used  to  describe  the 
explosive  and  was  employed  in  both  the 


theoretical  (DSD)  calculations  and  the  multi¬ 
material  simulations.  The  experiment,  theory  and 
simulation  were  found  to  be  in  excellent 
agreement  and  showed  that  for  a  large  class  of 
important  detonation  flows  one  can  use  the  DSD 
model.  DSD  assumes  that  the  detonation  shock 
propagates  along  its  normal  direction,  determined 
by  its  total  shock  curvature,  to  get  quantitatively 
accurate  results.  One  of  the  central  questions  that 
remains  is,  is  there  a  systematic  and  self- 
consistent  way  to  determine  the  location  of  the 
initial  detonation  shock  radius  (i.e.  the  location  to 
“light”  the  explosive)  to  start  the  DSD  simulation, 
as  a  precursor  to  a  program  bum  simulation? 

Therefore,  given  a  well-characterized 
initiation  source,  such  as  a  detonator  or  in  a 


computational  setting  a  prescribed  volume  of  high 
pressure  products  that  generates  a  shock  in 
unreacted  explosive,  one  needs  to  determine  the 
radius  for  the  starting  detonation  shock  shape  that 
is  consistent  with  the  assumptions  of  detonation 
shock  dynamics  (DSD)  theory.  The  same  study 
helps  determine  the  minimum  energy  and  size 
requirements  needed  to  achieve  a  stable  high- 
order  detonation.  We  also  are  able  to  assess 
significant  effects  due  to  shock  acceleration.  The 
inclusion  of  shock  acceleration  leads  to  one  of  the 
higher-order  theories  of  detonation  shock 
dynamics,  that  have  been  addressed  in  [4],  [5], 
and  [6]. 

In  what  follows  we  briefly  describe  a 
representative  explosive  application  that  was 
designed  using  a  PBXN-9  charge  with  19 
detonators  placed  in  an  equilateral  triangle  pattern 
array.  Figure  1  shows  a  framing  sequence  of  light 
captured  from  the  emergence  of  detonation  from 
the  initiation  points.  The  detonations  emanating 
from  the  initiation  points  do  not  initially  interact 
and  the  flow  can  be  modeled  as  spherical  or  axi- 
symmetric.  Later  the  initiating  detonation  shocks 
merge  and  the  flow  is  fully  3D. 

The  ideal  experiment  would  be  a  single  point 
initiation  characterization  into  a  PBXN-9  charge, 
but  the  19-point  experiment  was  a  related  and 
pertinent  data  set.  The  multi-point  experiment 
provides  redundancy  and  statistical  significance 
that  enables  an  estimate  of  the  initial  starting 
shock  radius. 


Figure  1.  Framing  sequence  a  19  point 
array  of  detonators  (order  proceeds 
from  top  left  to  bottom  right.) 

A  simple,  axi-symmetric  simulation  is  used  to 
estimate  the  effective  DSD  lighting  radius.  The 
lighting  radius  is  defined  to  be  the  minimum 
radius  from  an  initiation  site,  such  that  subsequent 
shock  evolution  according  to  the  shock  normal 


velocity,  total  shock  curvature  propagation  law  (a 
D-kappa  relation),  is  an  accurate  representation  of 
the  detonation  shock  motion.  The  properties  of 
the  initiation  system  also  determine  the  lighting 
radius  and  in  general  those  details  must  be 
considered. 

WIDE-RANGING  EQUATION  OF  STATE 
AND  RATE  LAW 

The  simulations  use  the  wide-ranging  equation  of 
state  and  a  reaction  rate  law  that  is  fully  described 
in  recent  papers  [3]  and  [1],  along  with  the 
calibration  procedures,  and  the  definition  of  the 
parameters  listed  below.  Basically  one  uses: 
Hugoniot  data  for  unreacted  explosive  and 
products,  the  assumption  of  a  simple  mixture  of 
reactants  and  products  and  some  information 
about  products  expansion  and  some  simple 
thermal  properties  to  calibrate  the  equation  of 
state.  Once  the  equation  of  state  is  fixed,  one  uses 
shock-to-initiation  data  (Pop-plot)  and 
experimentally  determines  D-kappa  data  for  the 
explosive  to  calibrate  the  rate  law.  The  resulting 
equation  of  state  and  rate  law  pair  is  then  used  to 
predict  behavior  over  a  wide  range  of  states 
without  further  calibration. 

PBXN-9  is  a  HMX  based  explosive  used  in 
conventional  munitions.  It  is  less  sensitive  to 
certain  insults,  but  is  still  a  high  performance 
explosive  (92%  HMX).  Unlike  PBX-9502,  there 
is  little  detonation  shock  velocity  versus  curvature 
data  available  for  PBXN-9,  with  the  exception  of 
overdriven  experiments  carried  out  by  L.  Hull, 
[7].  Calibration  of  PBXN-9  was  carried  out  in 
collaboration  with  B.  L.  Wescott. 

Products 

The  wide-ranging  equation  of  state  for  the 
products  has  six  (6)  adjustable  parameters.  The 
constraints  that  determine  the  parameters  are 
(quoting  from  ref.  [3]),  “...the  principal  isentrope 
must  pass  through  the  CJ-point  determined  from 
experiment.  Second,  the  principal  isentrope  and 
Rayleigh  line  must  be  tangent  at  the  CJ-state. 
Third,  the  total  energy  from  expansion  down  the 
isentrope  must  be  the  total  energy  liberated  by  the 
explosive.  Fourth,  adiabatic  gamma,  y,  at  large 
volume  must  be  the  ideal  gas  value  for  the 


products.  Fifth,  the  energy  from  the  expansion 
down  the  isentrope  must  be  partitioned  down  the 
high  pressure  part  above  0.1  GPa  where  useful 
work  -  driving  metal  is  done  and  the  low-pressure 
part  that  accounts  for  the  rest  of  the  chemical 
energy.  Sixth,  Gruneisen,  F,  must  be  adjusted  in 
the  high  pressure  region  to  describe  the  Hugoniot 
curve  for  the  overdriven  detonations  that  are  far 
from  the  principal  isentrope.  These  six 
requirements  determine  the  six  adjustable 
parameters.” 

This  calibration  process  showed  good 
agreement  with  the  JWL  EOS  expansion  data  of 
ref.  [8]. 

Reactants 

The  reactants  EOS  is  also  of  the  Mie-Gruneisen 
form  with  six  (6)  adjustable  parameters. 
Calibration  is  made  so  that  the  products  and 
reactants  Hugoniot  curves  do  no  cross  at  high 
pressures.  Experimental  shock  velocity  ( Us ) 
versus  particle  velocity  (Up)  data  of  ref.  [8]  is 
used  and  an  argument  about  the  reactant 
temperature  on  the  isentrope  is  used. 

A  mixture  equation  of  state  combines  the 
reactants  and  products  EOS  and  assumes  pressure 
and  temperature  equilibrium  conditions. 

Reaction  Rate  Laws 

L.  Hull,  [7]  carried  out  a  Mach  stem  experiment 
that  measured  the  D-kappa  relation  for 
converging  detonations  as Dn  =  8.56(1-0.7948  k) 
(k<  0),  where  Dn  is  measured  in  mm/ps,  and  k  is 
measured  in  mm"1.  At  this  time  no  experimentally 
determined  data  for  k  >  0  (divergent  geometry)  is 
available.  Therefore,  Hull’s  experimental 
correlation  from  the  Mach  stem  experiments  is 
used  to  constrain  the  theoretically  determined  D- 
kappa  relation.  Four  data  points  for  shock 
initiation  (PoP-plot  data)  available  from  ref.  [8] 
with  shock  input  pressures  in  the  range  from  3  to 
13  GPa,  were  used  to  constrain  the  rate  law. 

A  simple  depletion  law  with  a  pressure  dependent 
rate  was  assumed  of  the  form 


(  \ 

r  =  k  (l-A)v  (1) 

l  Pel) 

and  the  values  of  the  rate  law  parameters 
calibrated  with  Hull’s  data  and  the  shock 
initiation  data  are  shown  in  Table  1. 

TABLE  1.  Calibrated  value  of  parameters. 


k  (ns'1) 

N 

V 

14 

2.5 

0.7 

Once  the  wide-ranging  EOS,  and  rate  law 
description  is  specified,  one  can  readily  solve  the 
nonlinear  Eigen  value  problem  of  DSD  theory  to 
determine  the  normal  velocity,  total  shock 
curvature  relation  (i.e.  the  D-kappa  curve). 

INITIATION  OF  DETONATION  FROM  A 
CONSTANT  VOLUME  SOURCE 

To  test  the  transient  response  of  our  model  of 
PBXN-9  we  used  a  finite  volume  of  motionless 
products  of  radius  rinitiai  at  the  constant  initial 
volume  v0 ,  with  a  pressure  that  is  defined  by  a 
constant  volume  thermal  explosion.  This  is  the 
product  state  that  is  achieved  for  detonation  of 
infinite  velocity,  defined  by  the  intersection  of  the 
vertical  Rayleigh  line  with  the  product  Hugoniot, 
as  plotted  in  a  p,  v  -  plane.  The  initial  pressure  is 
approximately  15  GPa.  The  corresponding  run  to 
detonation  distance  for  a  similar  input  shock  to 
PBXN-9  is  approximately  2  mm.  With  this  initial 
state  well-defined,  we  only  varied  the  initial  (hot 
spot)  radius  rinitia^  and  then  solved  the  reactive 
Euler  equations  for  the  model  specified  above. 
The  lead  shock  location  is  recorded  along  with  its 
velocity  and  acceleration  as  a  function  of  time  for 
various  initial  (hot  spot)  sizes. 

The  results  of  such  simulations  are  shown  in 
Figure  2,  which  shows  the  lead  shock  normal 
velocity  Dn  as  a  function  of  the  lead  shock 
position  r.  Different  initial  (hot  spot)  radii,  that 
define  the  initiation  source,  are  explicitly  labeled 
for  r initiai  =  4.5  mm,  5.0  mm,  6.0  mm  and  7.0  mm. 
Note  that  the  D-kappa  relation  determined  from 
DSD  theory  is  also  plotted  as  a  D-r  relation, 
where  for  a  cylindrical  system  k  =  1/r.  Only  the 
top  branch  of  the  D-kappa  (D-r  relation)  is  shown, 
up  to  the  turning  point. 


Figure  2.  The  simulated  Dn-r  response 
for  PBXN-9,  initiated  from  various 
hot-spot  radii. 


Figure  3.  Normal  shock  acceleration 
(Dndot)  corresponding  to  Figure  2, 
versus  normal  shock  velocity. 


r-  [mm] 

Figure  4.  Normal  shock  acceleration 
(Dndot)  corresponding  to  Fig.  2,  versus 
the  shock  position. 


A  typical  transient  for  this  source  shows  a  slight 
initial  deceleration  followed  by  a  rapid 
acceleration  then  slower  acceleration  as  the  quasi¬ 
steady  detonation,  described  by  the  D-kappa 
relation  is  attained.  Figure  3  shows  a  plot  of  the 
shock  acceleration  versus  the  normal  detonation 
velocity  realized  for  various  initial  hot  spot  radii, 
as  seen  in  Figure  2.  Figure  4  shows  a  plot  of 
shock  acceleration  versus  shock  position  as  a 
function  of  initial  hot  spot  radii. 

The  results  show  clearly  that  the  D-kappa 
relation  ( D-r  relation)  is  an  attractor  of  the 
unsteady  solution.  In  these  simulations  one  can 
clearly  demark  a  radius  where  the  quasi-steady  D- 
kappa  response  is  achieved,  and  one  could 
consider  that  to  be  the  “lighting  radius”.  At  a 
sufficiently  large  radius  one  could  safely  assume 
that  the  detonation  propagated  according  to  the  D- 
kappa  rule.  But  it  is  clear  that  that  radius  depends 
strongly  on  the  strength  of  the  initiation  source 
and  the  subsequent  transient.  In  fact,  that  is  the 
main  point  of  our  article. 

Figure  5  shows  a  plot  of  the  “lighting  radius” 
as  a  function  of  the  “initial  starting  (hot  spot) 
radius”.  The  results  strongly  suggest  that  there  is 
a  minimum  radius  of  the  hot  spot  such  that  it  is 
not  possible  to  ignite  a  detonation  (below  4  mm). 


Initial  Starting  Radius  -  [mm] 

Figure  5.  The  lighting  radius,  (above 
which  a  DSD,  D-kappa  propagation 
rule  would  hold)  as  a  function  of  the 
initial  starting  radius. 


HIGHER  ORDER  DSD-THEORY  WITH 
SHOCK  ACCELERATION  EFFECTS 

For  an  ideal  equation  state,  Kasimov  and  Stewart, 
[6]  with  weak  curvature  asymptotics  and  slow 
evolution  on  the  time  scale  of  particle  transit 
through  the  reaction  zone,  derived  an  extended 
DSD  shock  evolution  equation  of  the  general 
form 

^  =  A(Dn)  ~  *  B(Dn) 

(2) 

where  A(Dn)  and  B(Dn)  are  functions  that  have 
explicit  forms  defined  by  the  EOS  and  rate  law. 
In  his  PhD.  thesis,  Kasimov  [9]  gave  a  similar 
development  for  general  equation  of  state.  Further 
Kasimov  and  Stewart  [6]  showed  that  for  states 
defined  by  blast  waves,  generated  by  explosives 
in  Hydrogen  Oxygen  mixtures  that  they  could 
predict  critical  energies  and  define  an  ignition 
seperatrix.  The  integral  curve  obtained  when 
plotted  in  a  D-r  plane  looked  in  a  qualitative 
manner,  quite  similar  to  those  shown  in  Fig.  2.  In 
the  absence  of  shock  acceleration,  the  D-kappa 
curve  is  obtain  directly  from  (2)  by  setting  the 
acceleration  term  equal  to  zero. 

The  general  theory  outlined  by  Kasimov  is 
being  applied  to  derive  the  result  corresponding  to 
equation  (2)  and  will  be  reported  on  shortly. 

CONCLUSIONS 

By  means  of  an  example,  as  applied  to  a  model  of 
PBXN-9,  prescribed  by  the  wide-ranging 
EOS/rate  law,  we  have  shown  that  significant 
transients  must  be  considered.  Further  the  lighting 
radius  for  the  application  of  DSD  must  take  into 
account  the  nature  of  the  initiation;  however  it 
seems  that  this  behavior  is  regular.  Future  work 
will  give  a  detailed  estimate  of  the  output  from 
detonators  and  companions  with  planned  new 
experiments  in  PBXN-9. 
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